Modeling the interaction of DNA with alternating fields 
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We study the influence of a THz field on thermal properties of DNA molecules. A Peyrard- 
Bishop-Dauxois model with the inclusion of a solvent interaction term is considered. The THz field 
is included as a sinusoidal driven force in the equation of motion. We show how under certain field 
and system parameters, melting transition and bubble formation are modified. 

PACS numbers: 87.14.gk, 87.50.U-, 87.16.A- 
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I. INTRODUCTION 

Terahertz (THz) technology and science have spread 
with ever-growing applications in military and security 
systems, medicine, biology, and researches. To exem- 
plify: security screening at airports [H, shipment in- 
spection 0, identification of concealed explosives, drugs 
and weapons 0, 0, cancer and burn diagnosis 
and in spectroscopy d, [H, [l^. Thus, knowing the ef- 
fects of THz radiation is critical for different scientific 
and technological purposes. 

Despite the presence of research on the biological ef- 
fects of THz radiation [ll|, there are still many contro- 
versies. In addition, the emergence of strong sources of 
THz radiation may contribute to the resolution of contro- 
versies over the mechanism of biological organization [l5l | . 
The potential of this perspective depends on the develop- 
ment of sophisticated pump-probe and multidimensional 
experimental techniques and the study of biological sys- 
tems in the controlled environments necessary for their 
maintenance and viability [l6| . 

In order to assess the possible genotoxicity of THz ra- 
diation with biological materials, in the framework of 
the THz-BRIGDE project [l3|, various studies were con- 
ducted to investigate the biological response. Recently, 
a careful analysis of these studies was successfully per- 
formed in [18], in which the authors explored the exis- 
tence of THz related effects on gene expression that can 
be unambiguously distinguished from thermal effects. It 
was suggested that THz radiation may affect gene ex- 
pression by perturbing the conformational dynamics of 
double-stranded DNA (dsDNA) [ll^. These studies 
were inspired by prior ones (22l . l23j|. As THz photons do 
not carry enough energy to directly alter chemical reac- 
tions, nonlinear resonance effects may cause local changes 
of breathing dynamics in these systems [24;, ^] . 

Motivated by this fact, Alexandrov et al. studied 
the influence of a THz field into the dynamics of a ho- 
mogeneous poly(A) DNA molecule with 64 base pairs 
(bps) [H, [20|. To model the interactions of dsDNA 
with THz field, they made use of the Peyrard-Bishop- 



Dauxois (PBD) model [2g] . They regarded periodic driv- 
ing and frictional terms in the absence of thermal noise. 
In that study, they found breather modes (localized peri- 
odic motions of the double strand) under certain condi- 
tions. Hence, they concluded that the main effect of THz 
radiation is to influence resonantly into the dynamical 
stability of the dsDNA. Though Swanson later showed 
that these breather modes can be eliminated by chang- 
ing the PBD model parameters or by including thermal 
noise [27j, he agreed that under assumptions concerning 
drag and drive forcing, breather modes can be generated 
at certain resonant frequencies. 

The PBD model is useful because it allows us to study 
the equilibrium and dynamic properties of DNA. It has 
potential ability to describe the melting transition and 
denaturation bubbles of the dsDNA such as those that 
occur during the initial stage of the transcription pro- 
cess. Some specific correlation with transcription ini- 
tiation sites has been claimed [12, [H-[3l| and debated 

mm. 

Tapia et al. [sj] set suitable parameter values for PBD 
model and studied the formation and stabilities of bub- 
bles in the system. They used an enhanced model that 
includes solvent interactions through the addition of a 
Gaussian barrier to the Morse potential [s^. This bar- 
rier modifies the melting transition and the dynamics of 
the molecule. They also focused on the application of the 
Principal Component Analysis (PCA) of the trajectories 
under equilibrium conditions. 

Even if both the PBD model and the interaction with 
external field taken as Alexandrov are quite simple, they 
could give some insights in the understanding of the in- 
teraction of THz field with the DNA molecule. For these 
reasons our purpose is to study the effect of the THz 
field on thermal and dynamic properties of DNA: melting 
transition and bubble formation at physiological temper- 
atures in the framework of the PBD model. We used a 
modified version of the model with the inclusion of sol- 
vation barrier and included thermal noise. Both homo- 
geneous and heterogeneous chain are studied in order to 
approach more to reality. The heterogeneous chain is the 
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adeno-associated viral P5 promoter (AAV P5 promoter 
or P5 promoter), which has been widely studied [22l. [29| 
and plays an important role in the AAV DNA replication 
[36'] and the regulation of the AAV gene expression [STj . 

This work is divided as follows: in section II we de- 
scribe the model; in section III the methods are summa- 
rized; in section IV an analysis of the response of the sys- 
tem at different frequencies and field amplitudes is made. 
After that, we will show the influence of THz field, with 
specific parameters, into melting transition (section V) 
and denaturation bubbles (section VI). The last section 
is devoted to conclusions. 



We used the same value parameters of the PBD model 
with solvation barrier that appear in reference [s^ : 
Dat = 0.05185ey, Gat = 0.1556el/ , yoAT = O.sA, 
bAT = 0.03125l^ K = 0.03eFl^ p = 3 and ^ = 
0.8l"\ 



III. METHODS 

In order to study the behavior of the system we have 
performed molecular-dynamics numerical simulations of 
the Langevin equation: 



II. THE MODEL 

The PBD model is a mesoscopic dynamical model of 
the DNA molecule. It describes the stretching of the 
bonds between the bps through a single variable, which 
condenses all the atomic coordinates of a bp. This model 
ignores the helicoidal structure and uses the Morse po- 
tential to model hydrogen bonding between bps. A non- 
linear inter-pair stacking potential is also considered. We 
use a modification of the PBD model including a solva- 
tion barrier inside of the Morse potential. This barrier 
prevents the closing of the base once it is opened. The 
total energy of the system is then approached by: 



2m 



+ Viyn) + W{yn,yn-i) 



(1) 



In this equation V{yn) is an on-site potential that de- 
scribes the interaction between the two bases of a pair. 
It is represented by the Morse potential and a Gaussian 
barrier is added: 



v(y) = Die-^y - 1)2 + Gc^'^y-y"^'/''. 



(2) 



Parameters D, a, G, yo and b are sequence dependent. 
Following [ssj, in our simulations we will use Dec = 
1.51? AT and acG — l-5aAT- D is the bp dissociation 
energy and a sets the amplitude of the potential well. 
The barrier height is controlled by G, its position and its 
width are given by yo and b respectively. A reasonable 
selection for such parameters is G = 3D, yo = 2/ a and 
b = \l2o? [3^. The term W{yn,yn-\) accounts for the 
stacking interactions and is given by 



m^^ + m — = 9[W{yn,yn+i + W(yn-i,yn)] 
df^ dt dyn 



dV 

dyn 



+ £,n{t)+ Acos(Lot), (4) 



where m is the mass of the bp, 7 is the effective damp- 
ing of the system and ^(i) accounts for thermal noise, 
(6.(0) = and (en(i)Cfc(i')) = 2m^kBTSnkSit - f), 
T is the bath temperature, A and oj are the field 
amplitude and frequency respectively. The equations are 
numerically integrated using a stochastic Runge-Kutta 
algorithm [4l|, |42| . We thermalize during 200ps without 
field and 200ps with the field before any record is 
made. Simulations of melting transition are performed 
using periodic boundary conditions while those of 
thermal bubbles with fixed boundary conditions as 
in reference [33 |. The P5 promoter is given by the 
69 bps: "5- GTGCCCATTTAGGGTATATATGGCC- 
GAGTGAGCGAGCAGGATCTCCATTTTGACCG- 
CAAATTTGAACG - 3" . For methodological issues we 
analyze homogeneous chains as well. In these cases the 
chains have the same number of bps as the P5 chain. 

To show the influence of field parameters in the system 
response we use the mean displacement (y) defined as 



(5) 



n,t 



where N is the number of bps and ts is the simulation 
time. For the melting transition we also measure the 
mean energy (u) as a function of the temperature 



W{yn,yn-{) = iif(l + pe-^(^"+^'-))(y„-2/„_i)2. (3) 

The effect of this term, whose intensity is governed 
by p, is to change the effective coupling constant from 
K{l-\- p) to K when one of the bps is displaced away from 
its equilibrium position. The parameter 5 sets the scale 
length for this behavior. Alternatively, inhomogeneous 
stacking energy can be also considered [39|. 



In order to calculate the opening probability and life- 
time of a bubble we followed reference [l^. The prob- 
ability Pn{l,tr) for the existence of a bubble of certain 
length / of bps, threshold tr and beginning at n*'* bp is 
calculated as: 
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Pnil,tr) 



- T. 



At[q^{l,tr)] 



(7) 



where qn{l,tr) counts for the bubbles of duration 
At[qn{l,tr)]. The average bubble duration r is calcu- 
lated as the average time of a bubble of a given shape 
over all occurrences of that bubble, 



At[qn(l,tr)] 



[<ln{l,tr)] 



(8) 



We can extract information from a large set of data in a 
multidimensional phase space through the PCA. It allows 
to reduce the dimensionality of the variable to those that 
include most of the fluctuations of the original system 
[40| . From an operational point of view, we have to build 
the N X N correlation matrix. So, 



(9) 



The diagonalization of this matrix allows us to obtain 
an ordered set of eigenvalues (Ai > A2 > A3...) with their 
corresponding eigenvectors (wi, W2, W3....). The amount of 
fluctuations is given by the eigenvalues. The new coordi- 
nates are ordered in such a way that most of the system 
fluctuations are retained by the few first ones. 



IV. SYSTEM RESPONSE FOR DIFFERENT 
FREQUENCIES AND FIELD AMPLITUDES 

Before performing the simulations for the melting tran- 
sition and bubble formation, some preliminary steps 
should be done. First, we look for the frequency val- 
ues at which maximum responses are obtained for each 
sequence. These values depend on temperature. At low 
temperatures they should be in the order of the linear 
resonances of the system. Nonlinear oscillations become 
important at intermediate temperatures while above the 
melting temperature, frequency values belong to those of 
a Gaussian chain. 

The system dynamics depends on damping values as 
well. Thus, we have considered two values for the damp- 
ing coefhcient 7 = Ips^^ and 7 = 9.8ps^^. Figure [1] 
shows the behavior of three sequences with 69 bps. The 
calculation has been performed with a field amplitude 
A = 50pN, damping factor 7 = lps~^ and two tempera- 
ture values: T = 210K and T = 290K. 

Maximum responses occur at certain frequency values, 
even for large field amplitudes (not shown) . Some modes 
are activated when temperature increases, (see figH]). 
The frequency bandwidth slightly increases when tem- 
perature rises. These bands are around w = 9rad/ps and 
oj = Vlrad/ps for AT and CG chains respectively. The 
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FIG. 1: (Color online) Frequency dependence of {y) with A = 
50pN and 7 = lps~^.a) AT chain is represented by lines with 
triangles and CG one by circles. Filled markers are for T = 
290K and empty are for T = 210^". b) P5 promoter. Solid 
line is for T = 2107^ and dashed line for T = 290K. Figure 
inset: Frequency dependence of {y) for each base at T = 
2907^. 



increased of (y) for certain frequency values could lead 
to bubbles formation or the full melting. The resonant 
frequencies of P5 promoter are the same as those for the 
AT and CG chains. In other words, P5 promoter behaves 
as if it would be composed by two homogeneous chains 
on the frequency space. This behavior may be under- 
stood because AT and CG bps number on the chain are 
approximately equal. However, the (y) values are less 
than the corresponding to the homogenous sequences. 
At T = 210K the maximum response for P5 promoter 
occurs around uj = 9rad/ps. Only the AT bases are 
stimulated and for this reason localized openings are ob- 
served around AT richer regions, (see the inset of the fig. 
[T]b). At T = 290K the maximum response occurs at 
both io — 9rad/ps and uj — llrad/ps. The CG bases 
are also stimulated. In this case, there is enough energy 
for opening AT bases as well due to stacking interactions 
and the whole chain opens. Figure [5] illustrates the fre- 
quency dependence of {y) with 7 — 9.%ps~^. With a 
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FIG. 2: (Color online) Frequency dependence of {y) with ^ = 
lAApN and 7 = 9.8ps~^. a) AT chain is represented by lines 
with triangles and CG one by circles. Filled markers are for 
T = 2<dQK and empty are for T = 21QK. b) P5 promoter. 
Solid line is for T = 21QK and dashed line for T = 290A'. 
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FIG. 4: (Color online) Amplitude dependence of (y) with 
7 = 9.8ps~^. a) AT chain is represented with triangles and 
CG one with circles. Filled markers are for T = 290/s' and 
empty are for T = 21QK. b) P5 promoter. Squares are for 
T = 210K and circles for T = 2'dQK. Lines are guide to eyes. 
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FIG. 3: (Color online) Amplitude dependence of (j/) with 
7 = lps~^ . a) AT chain is represented with triangles and 
CG one with circles. Filled markers are for T = 290A' and 
empty are for T = 210A'. b) P5 promoter. Squares are for 
T = 210-?^ and circles for T — 290K. Lines are guide to eyes. 



larger damping, it is no longer possible to obtain res- 
onant frequencies because the stochastic term becomes 
dominant. 

We have also calculated the amplitude dependence of 
(y) to determine A values for the next simulations. The 
frequency values of maximum responses are used: uj = 
9rad/ps for AT and P5 chains and uj = V! rad/ps for CG 
chain. Results for cj — llrad/ps for the P5 chain are 
similar to those oi uj — 9rad/ps. Figures [3] and |4] show 
the results for 7 = Ips^^ and 7 — 9.8ps~^ respectively. 

According to fig. [3] for 7 = lps~^, we can use A = 
10, 25 and 50pN for the three chains. In the case of 
7 = 9.8ps~^ the values are A = 50, 144 and 200pN. 



V. MELTING TRANSITION 

This section focuses on the study of the melting tran- 
sition of the homogeneous chains and P5 promoter. Our 
goal is to analyze how a THz field modifies melting tem- 
perature Tjyi and the transition width AT. These param- 
eters were determined in reference (33 | for the uniform 
chain of AT bps without external field. Following the 
same criteria we can determine them when the external 
field is applied. We calculate the mean potential energy 
and mean displacement as a function of the tempera- 
ture. Let us determinate two temperatures as follow: the 
temperature Ti estimates the beginning of the transition 
where {y{T)) crosses 0.5 A. The larger one, T2, provides 
the onset of the linear behavior in {y{T)). Both quanti- 
ties are defined in the following form: Tm = (Ji -I- T2)/2 
and AT = (T2 — Ti)/2. Figure [5] displays melting transi- 
tion curves for the AT chain and P5 promoter. Melting 
transitions curves for CG and P5 chains at w = Ylrad/ps 
have similar behaviors (not shown). 

While the mean potential energy curves coincide above 
T2 for different field amplitude and damping values, the 
mean displacement curves differ. This can be explained 
because the action of the field is such that the difference 
{y-n — Un-i) between two successive bps at each time step 
remains constant. In this region the potential energy 
depends only on the difference (j/„ — yn-i)- Figure [5] 
shows the behavior of the melting transition temperature 
versus field amplitude for the three chains. 

The differences between the melting transition with 
and without applied field are remarkable (see figs. [5] and 
0. Due to the applied external field, the chains melt at 
lower temperatures. This fact has been already noted by 
Swanson [27| for a homogeneous chain. Differently from 
this work, the values of Tm we obtained here are larger 
than the one there reported. Behavior of the transition 
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FIG. 5: (Color online) Melting transition at the frequency 
cj — 9rad/ps. a) AT chain with damping factor 7 = lps~^ , 
b) AT chain with damping factor 7 = 9.8ps~^ , c) P5 chain 
with damping factor 7 = lps~^ and d) P5 chain with damping 
factor 7 — 9.8ps~^. 
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FIG. 6: Melting transition temperature, a) damping factor 
7 — lps~^ . frequency uj = 9rad/ps for AT and P5 chains 
{P5ujAt) and cj = llrad/ps for CG and P5 chains (P5[^cg)-, 
b) damping factor 7 = 9.8ps~^. frequency uj = 9rad/ps for 
AT and P5 chains and cj = llrad/ps for CG chain. 



width is more complex because it depends on the num- 
ber of bps that has been opened at certain temperature. 
The field allows both opening and closing of bps. For 
7 = 9.8ps~^ the effect of the field decreases and larger 
values of field amplitude are needed to lower the transi- 
tion temperature. 



VI. BUBBLES FORMATION 

Next, we study the THz field infiuence in the bub- 
bles formation at T = 2%QK with the parameters cho- 
sen in previous section. This temperature belongs to the 
premelting range, in which it has been shown that the 
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FIG. 7: (Color online) Probability opening distribution P„ 
with frequency cj — 9rad/ps, damping factor 7 — Ips^^ and 
T = 290^". a) AT chain b) P5 promoter. 



highest opening probability coincides with important bi- 
ological sites like the start transition site (TSS) and the 
TATA box (28l - [30| . In order to avoid unphysical denatu- 
ration process due to finite-size effects, we add 10 CG bps 
sequence to the ends of the P5 promoter to create hard 
boundaries. The extremes are set to zero , avoiding the 

Figures [7] and [8] show 



complete opening of the chain 
the bubble size and lifetime distribution respectively for 
the AT chain and P5 promoter with uj — 9rad/ps and 
7 = lps~^. Pn and T„, given by equations [7| and [51 are 
defined with a threshold value of tr = l.bA. These mag- 
nitudes are represented as a function of the bubble length 
and index site. The opening probability and bubbles life- 
time are given in color scale. In these figures the 10 bps 
at the beginning and the end of the sequence are not in- 
cluded. The -1-1 in Base Pair Index refers to TSS position 
in P5 promoter. In the homogeneous chain there is no 
TSS but we keep the same notation for convenience. 

Figure [7] shows how the external field enhances the 
opening probability at the frequency value u = 9rad/ps. 
Results agree with those of the melting transition. With- 
out field, the larger probabilities in P5 promoter occur in 
two sites of biological interest as previously reported: the 
TSS (represented by -f 1) and the TATA box (between - 
30 and -40). Increasing the field amplitude leads to the 
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increase of opening probability at these sites and helps 
the opening of others. Without the THz field, the open- 
ing probability of P5 promoter at the TATA box is higher 
than at the TSS [Ulil]. The most persistent bubbles 
are found at the sites that have been pointed out before, 
see Fig. [8l In our simulations, the lifetime values depend 
on the selection of the parameters of the modified PBD 
model. Thus, the results in some cases can be different 
respect to those reported in literature when the THz field 
is not applied. For the AT chain, the field makes bub- 
bles more stable but for the P5 promoter it does not. In 
this case, the heterogeneity of the chain plays a crucial 
role, because we have used to — 9rad/ps. The decrease 
of bubbles lifetime may be explained because the energy 
of the field favors both opening and closing events of the 
bps. The opening and closing kinetics is governed by 
the solvation barrier and the applied field. With bar- 
rier and without applied field, the kinetics is controlled 
by the presence of two equilibrium states separated by 
the solvation barrier. Closing events are more difficult 
and bubbles live longer than those in which the barrier is 
not considered jsj]. In the other hand, at this frequency 
CG bps are not stimulated. If we use uj = 17rad/ps, 
CG bps are stimulated and the whole chain opens, as ex- 
plained previously. Fig. |9l The same calculations were 
performed with 7 = 9.8ps^^ (not shown). As in previous 
results, damping factor modifies the chain dynamics and 
high field amplitude values are needed to open the chain. 

Finally, we compare the results obtained from the PC A 
with the average displacement values of each base (see 
figure [T0| . The 10 CG bps at the ends of the sequence 
are included in the figures. We only show the results for 
P5 promoter with a field amplitude oi A = 50pN. 

A good correlation between top figures and bottom 
ones is verified. Localized eigenvectors span over regions 
of nine bps, which fairly correspond to the width of the 
bubbles, as in reference [33|. For other field amplitudes 
and the homogeneous chains there is a good correlation 
too. 

Our results suggest that the THz field influence can be 
viewed in two main directions. First, the applied field 
tends to facilitate melting and bubble formation which 
could in principle affects processes like transcription or 
replication. The driving forces needed for observing such 
effects are large compared with those found in physical 
reliable conditions for in vivo exposure [27]. We do not 
disagree with this conclusion. We have just used a simple 
model for such a description and the magnitudes values 
may no match the real ones. On the other hand, the 
external field could be used to detect biologically signifi- 
cant sites by increasing the opening probability at these 
sites without the melting of the chain. In this scenario 
the effective drive could be larger. 
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FIG. 8: (Color online) Average lifetime distribution Tn with 
frequency u) = 9rad/ps, damping factor 7 — lps~^ and T = 
2mK. a) AT chain b) P5 promoter. 
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FIG. 9: (Color online) Probability opening distribution P„ 
and average lifetime distribution r,i for P5 promoter with fre- 
quency cj = VJrad/ps, damping factor 7 — lps~^ and T = 
2907^. Top figure: A = 25pN. Bottom figure: A = 50pN. 



VII. CONCLUDING REMARKS 

We have studied in the framework of PBD model the 
infiuence of THz field on homogeneous chains and the 
heterogeneous P5 promoter. Thermal properties of these 
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FIG. 10: (Color online) Top figure: probability of opening for 
P5 promoter. Bottom figure: the three first PC eigenvectors 
corresponding to eigenvalues Ai, A2 and A3 at T = 290K. a) 
7 — lps~^ . b) 7 = 9.8ps~^. 



sequences have been studied by including thermal noise 
and a solvation barrier in the model. The influence 
of THz field depends strongly on field parameters (fre- 
quency, field amplitude) and system parameters (poten- 
tial parameters, damping). In spite of previous results, 
we do not obtain breathers modes (oscillatory solutions) 
but rather we find that ac-fields favors the formation of 
bubbles. 



We have also identified the frequency resonant bands 
that mostly increase the opening of the chain. The posi- 
tion of the bands are sequence dependent and distinguish 
the AT-rich from CG-rich regions. This could increase 
experimental resolution in order to detect sites like TSS 
if small field amplitudes are used. More study in this 
direction needs to be done, for instance considering more 
complex interaction between bases and the external field. 



Finally, we have numerically obtained that the PCA 
can be also used to get information for out equilibrium 
systems. 
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